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Abstract 

ATUS-PRO is a solver-package written in C-l—I- designed for the calculation of numerical solutions of the 
stationary- and the time dependent Gross-Pitaevskii equation for local two-particle contact interaction utilising 
finite element methods. These are implemented by means of the deal.II library mm- The code can be used 
in order to perform simulations of Bose-Einstein condensates in gravito-optical surface traps, isotropic and full 
anisotropic harmonic traps, as well as for arbitrary trap geometries. A special feature of this package is the 
possibility to calculate non-ground state solutions (topological modes, excited states) [3l|4l[5] for an arbitrarily 
high non-linearity term. The solver-package is designed to run on parallel distributed machines and can be 
applied to problems in one, two, or three spatial dimensions with axial symmetry or in Cartesian coordinates. 
The time dependent Gross-Pitaevskii equation is solved by means of the fully implicit Crank-Nicolson method, 
whereas stationary states are obtained with a modified version based on our own constrained Newton method 
[5]. The latter method enables to find the excited state solutions. 
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1 Introduction 

In recent years, since the first experimental realisation in 1995, Bose-Einstein condensates (BEG) have caught 
increased attention in the experimental as well as in the theoretical physics community. The properties of BECs 
can be well described by means of a partial differential equation of the type “nonlinear Schrodinger equation” 
(NLSE). Here this is the Gross-Pitaevskii equation (GPE). In this description, the BEG is treated as a non- 
uniform, interacting Bose gas at zero temperature. The term “interacting” refers in the GPE-description to at least 
two-particle interactions, which are modelled in a mean-field approximation and give rise to a non-linear term. 

There exist only a few analytic solutions of the GPE, like for example soliton solutions in one dimension, the 
frequently used Thomas-Fermi approximation for ground states and variational approaches [^. Also on basis of 
the Thomas-Fermi approximation the dynamics of BEGs in time-dependent traps can be described by means of 
the scaling approach [7]. However, the approximations and assumptions made in these approaches are not strictly 
fulfilled in certain situations. Especially when it comes to high-precision measurements, deviations from idealised 
setups have to be taken into account. For example, the Universality of Free Fall (UFF) is to be tested with quantum 
matter [Hiiiiin] where possible violation signals would be of extremely low magnitude. Since only by including all 
relevant error sources in the calculations, one is able to distinguish these signals from systematic effects introduced 
by the experimental equipment or, in general, by the environment. For this, three-dimensional simulations have to 
be performed that naturally lead to a high demand of computational power. Often, one needs to utilise massively 
parallel computer systems. 

For example, a comprehensive program package for the simulation of the Gross-Pitaevskii equation has been 
developed by Muruganandam et al. m and its extension to G by Vudragovic et al. [HITS]. It allows to solve the 
stationary and non-stationary case for ID, 2D and 3D for fully anisotropic traps. 

In our recent investigation [3] we developed a contrained Newton method in G-l—h by which numerical solutions 
for excited states in gravitational surface traps and harmonic traps can be obtained. This was implemented for 
one-dimensional stationary problems using finite differences. 

Goncerning physical motivation to study excited states it has to be mentioned that Biicker and coworkers 
demonstrated experimentally the vibrational state inversion of a Bose-Einstein condensate nans]. The experiments 
were guided by simulations based on optimal control theory which is a dynamical method. There, the toolbox OGT- 
BEG |16] was used in order to get experimental prescriptions how to control the time-dependent trapping potential. 
By means of this procedure, the first excited state could be reached. As a test of how particular higher order excited 
states presented in [3] could be realised experimentally, we performed some preliminary calculations utilising the 
OGT-BEG package. We were able to specify time dependent external potential configurations from which we 
obtained the excited states up to the third mode. This shows that the predicted states could be produced, in 
principal, experimentally. The results will be published elsewhere. 

Excited states also find applications in the frame of exploring the interaction of quantum matter with gravity. 
Energy eigenstates of ultracold neutrons in a gravitational trap have been investigated in order to test Newton‘s 
inverse-square law at small distances by Abele et al. m- In their setup neutrons fall in the linear gravitational 
potential and bounce off a neutron mirror. The neutron flux is measured as a function of the absorber height above 
the mirror. Therefore, it is a function of the vertical discrete Energy component En and thus depends on the mode 
number n. By means of this, a constraint in the range between I ^m and 10 /im could be obtained for Yukawa-like 
potentials. As a side remark, such setups are also known as “atom trampolines” or “quantum bouncers” whose 
solutions for the linear case, i.e. the Schrodinger equation, are given by Airy functions |18| . Furthermore, it has 
been suggested to use quantum bouncers in order to complement tests of the Universality of free fall , see m 
and references therein. The main conclusion is that experiments which investigate the dynamics of wave packets, 
independently of the complexity of the initial state, always probe the ratio mg/rrii of the gravitational rUg and 
inertial mass rrii. However, the probability density can depend on the product of both masses like or 

other combinations. Moreover, the discrete energy spectrum can behave according to the scaling law . 

This could be explored, in principle, by means of spectroscopic methods. This would add to existing experimental 
procedures and setups for comparing rUg and m^. 

Here we present a highly improved version of our previously mentioned G-l—h code which has been extended to 
two- and three-dimensional setups and comes with the following features: (i.) calculation of ground state solutions, 
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as well as (ii.) stationary excited states, (iii.) dynamics simulations via real-time propagation, and all calculations 
are based on (iv.) finite element methods (FEM). By the latter, the code can now be used as a basis for solving 
problems of higher spatial complexity on massively parallel computer systems (of course also on single multi-core 
systems), e.g. time propagation of highly excited states in anisotropic time-dependent anharmonic traps. 


2 Physical model and notes on the algorithm 

2.1 The model 

The physical properties of Bose-Einstein condensates are well characterised by means of the Gross-Pitaevskii equa¬ 
tion which in its time-dependent and dimensionless form can be written as 

idt^{x,t) = (-A + V{x) + 7 |\E'(x,f)|^) (1) 

It describes the dynamics of a many body bosonic system with local self interaction at zero temperature. Eor 
this equation to be a good approximation, it is assumed that the s-wave scattering length is much less than the 
mean interparticle spacing. The wave function ^(x,t) is complex and |l'(xjt)p is interpreted as the local density. 
The local self interaction strength is given by 7 £ IR (which depends on the s-wave scattering length) and its sign 
(positive or negative) determines whether the particles of the condensate repel or attract each other. Usually, the 
external potential V(x) models a trap in order to spatially confine the condensate, but can also account for e.g. 
external perturbations on the system. Here we assume that it is bounded from below, i.e. V(x) > 0. 

By means of the following ansatz 


(x, t) = (j){x) exp (i/rt), 4>{x) € R, (2) 

space and time are separated and one obtains the stationary Gross-Pitaevskii equation: 

(-A + V(x) - fi + j(/>'^(x)) (p(x) =0, (3) 

where /r is the dimensionless chemical potential. Note that all stationary solutions are real. 

As default setting of the atus-pro package, a gravito-optical surface trap (GOST) is utilised in two dimensions. 
That is, the gravitational attraction is superimposed by a harmonic trap, i.e. 

V(x,p') = ^x+ uj^y^, (4) 

for Gartesian coordinates and 

V(r,z)=uj^r^ + /3z, (5) 

for the axially symmetrical case, where r = Eor the three-dimensional case the potential of the GOST 

is implemented for Cartesian coordinates according to 

V{x,y,z) = Px + ixly'^ +ujIz‘^. ( 6 ) 

Inserting an infinite repulsive potential wall at a; = 0 or z = 0 establishes the boundary condition at the bottom of 
the trap, whereas (j){r = R) = 0. This ensures that the Hamiltonian is bounded from below. Then the potentials 
([5|1 and ([S]) model the behaviour of a bouncing quantum system, i.e. the previously mentioned “atom trampolines” 
or “quantum bouncers”. In our package fully anisotropic harmonic traps are also predefined, that is 

V{x,y) =ujIx^+ ujly'^, (7) 

and 

V{x,y,z) =ujIx‘^ +ujly'^ +ujIz'^. ( 8 ) 

The dimensionless parameters P and uj,u!x,ujy,ujz represent the gravitational acceleration and the angular frequen¬ 
cies of the harmonic trap. 

2.2 The contrained Newton method for the stationary case 

Solutions to the stationary Gross-Pitaevskii equation are obtained by means of a constrained Newton algorithm 
which is a slightly modified version of the algorithm presented in [5]. In our original implementation which was 
used to obtain the results in [3], we used the Finite Difference Method (FDM) for spatial discretisation. However, 
in order to release a public version of this code we decided to reimplement the algorithm using more sophisticated 
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Figure 1: Enhanced Newton method: Flow chart of the outer loop, s is outer loop counter, k is the inner loop 
counter. 


methods. In particular, we make use of libraries for Finite Element Methods (EEM). The flow charts in figure [T] 
and figure [H explain how the stationary solutions by means of the breed programs are calculated. 

We start reformulating the problem by means of multiplying the stationary GPE with a test function 
h{x) G and integrating over the volume 

f [—A + V{x) — ^ + (j){x)h{x)dV = 0, (9) 

Ja 

which by means of Green's identity becomes 


(V(/>(£) • Vh{x) + {V{x) — /i) (j){x)h{x) + 'y(f>^{x)h{x)) dV = 0. 


( 10 ) 


Note that terms at the boundary did vanish because of the fact that h{x)\gQ is zero. The expression (uni) is also 
known as the weak formulation of the partial differential equation @. Solutions 4>{x) are critical points of the 
corresponding functional 


A-W^lA = + ^{V{x) - p)(l)'^{x)+ ^(l)‘^{x)^ dV 


( 11 ) 


and determined in our algorithm by solving equation (fTUl) instead of the stationary GPE ([3]) . That means we look 
for the zeroes of cni) which are not unique. In fact, for ^ oo one gets infinitely many zeros m US [11121111]. 
However, for a fixed /i there exist only a finite number of zeroes. In contrast, if we fix the particle number then 
we have a discrete spectrum of infinitely many different solution. The solution with the smallest p, is then called 
“ground state” of the system and belongs to a minimum. All other solutions are min-max saddle points, like in the 
linear case, too. The excited states, or non-ground state solutions, are of that kind. 

Before starting the Newton iteration loop, the initial wave function is assembled as a linear combination of two 
functions. 


is,k 

Kef 


= t: 


s,k 


>o + tl’K 


s.k 


( 12 ) 


where in the first step the coefficients and tf^ are set automatically or by means of a determined value in 

the parameter file params. prm. If set to automatic initialisation, tQ^~^ = tf^~^ = (for a better overview. 
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please see the flow chart in figure[T]). Here and in the following s and k denote iteration counters. Note that initially 
the function is set to zero. In contrast, (/>g is chosen either (i.) automatically w.r.t. the desired nodal 

structure of the solution, or (ii.) the function can be specified explicitly in the parameter file. In the first case, (/)g 
is given by the solution of the corresponding linear problem, i.e. the solution of the Schrdinger equation with the 
same potential. 

For example, in the case of a three dimensional harmonic trap the potential is completely characterised by 
means of the frequencies (omega_x,omega_y,omega_z) given in params.prm. In case (i.) (j)^ is then determined by 
specifying the quantum numbers (QN1_x, QNl_y, QN1_z) in the same parameter file, which refer to the stationary 
solution of the Schrddinger equation. By means of this, excited state solutions can also be calculated. The precise 
form of the implemented eigenfunctions will not be reproduced here, in fact this can be looked up in the doxygen 
documentation which is included in the atus-pro program package. 

Concerning case (ii.) all functions and combinations thereof which are defined in the C standard math library 
are allowed. However, the manually determined (/)g has to respect the boundary conditions. 

It is important to mention that after the first Newton iteration step the coefficients and are fixed by the 
constraints (ITT)) and (IT^ . 

The FEM-form of equation m reads 



The integral over the volume translates into a double sum over all cells K and all quadrature points in each cell 
times a weight , given by the product between the Jacobian of the cell and the weight of the quadrature formula. 
Here the indices are: i denotes the degrees of freedom (in FEM language these are the expansion coefficients in 
front of the shape functions), K is the sum over all cells and q signifies the sum over all quadrature points in each 
cell. hi{x^) is the shape function belonging to the i-th degree of freedom. 

From this the Jacobian is derived by means of taking the first functional derivative, hence 



This is needed in order to calculate the search direction of the algorithm for the Newton method, namely 


jS, k iS, k 77Tc 

Jy d, =F, 


s,k 


must be solved and next the function 


+1 = - sgn{tf)F^\ 


(15) 

(16) 


is calculated. 

It is important to note that the sign of the coefficient is taken into account which is due to the fact that the 
problem has Z 2 symmetry. Since a solution (/) of the Gross-Pitaevskii equation can be replaced by —(f) being also a 
solution, this leads to a change of sign for the coefficients as well, i.e. —tg ^ and —In this case it happens that 
—>■ —however maintains its sign. This would induce a sign change of d®’^ which is compensated by 
sgn(ti’''). 

In the next step of the Newton loop the constraint of our Newton algorithm is specified by the following system 
of equations 


V(j)lf/\x^)V(j)^o{x^) + {V- ) +7 )) 

K,q L 


K,q 


= 0, (17) 
= 0. (18) 


The mathematical meaning of both equations is to obtain a special point (j^ref function space where the 

F 2 -gradient given by the l.h.s. of time independent Gross-Pitaevski equation (l3|) vanishes in the direction of (j)o 
4>i^. That is, one follows all these points until a critical point is found at which all derivatives vanish. 

Upon finding new values and from the k-th. step, a new reference function is assembled by means 

of equation ca) and the iteration loop starts again until the L 2 -norm of the F 2 gradient reaches a predetermined 
value. That is, the F 2 norm of the left hand side of equation © should be close to zero. 
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Figure 2: Enhanced Newton method: Flow chart of the inner loop, s is the outer loop counter, k is the inner loop 
counter. 


A typical simulation run can now be performed in two ways: (i.) the code will compute a large number of 
solutions - each for a different chemical potential /i® in increasing order or (ii.) one single solution for a fixed /i is 
calculated. 

In case of (i.) the code will try to find a solution for an initial chemical potential and upon successful 
completion of the Newton iteration loop it proceeds to determining the next solution for a new = /r® + A/r. 
That is, the Newton loop is evoked each time a new /i® is fixed and it runs as long as s < Ndmu. Here, the step size 
A/x is specified in the parameter file params. prm. We have to emphasise that this way of calculating solutions is 
well suited for computations similar to those in [3]. 

As a remark, solutions are currently restricted to real functions. Moreover, the reader should be aware of 
the fact that the solutions are not normalised to one. Due to the one-to-one correspondence between fj, and 7 (see 
for example Fig.8 in m) the final particle number N is determined by this pair. In order to get solutions normalised 
to one, this has to be done by hand and the nonlinearity 7 must be adjusted accordingly, i.e 7 —> 7 = A7. This 
will be addressed in subsection 13.61 

Here, a general statement is at hand with respect to how degenerated states solutions are handled. Excited 
state solutions are degenerated by rotational symmetry in case of isotropical harmonic trapping. Usually, this would 
lead to a failure of the standard Newton algorithm. However, since we use a special constraint, this problem can 
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be circumvented and the Jacobian of our Newton algorithm is not degenerated. This translates to the fact that 
the structure and the positions of the nodes are fixed by the choice of </>[]. Thus, once chosen, in the case of a 
two-dimensional harmonic trap, the positions of the lines of nodes in the x — y-plane will not change in the course 
of the calculations (see figure|6]). In order to obtain rotated solutions w.r.t. the previous one, one has to rotate <j)Q, 
too. 

However if /x is accidentally chosen in such a way that it corresponds to a bifurcation point, then the Jacobian 
becomes degenerated. In this situation the program terminates with a note that the matrix is singular. 

For a further and more detailed description of the algorithm, we recommend you to consult the doxy gen docu¬ 
mentation. 


2.3 Notes on the algorithm for the time-dependent solutions 

In order to advance the complex wave function 'I' from time t to t + At we have implemented the fully implicit 
Crank-Nicolson method for the time dependent Gross-Pitaevskii equation. This scheme is unitary in time, conserves 
the total number of particles and the total energy. 

In order to derive the system of equations, we start with the following ansatz 


’'At I 


V- 


T(x, t -f At) + 'i'(x,t) 


/ 'I'(x,t + At) + 'i'(x,t) 

V 2 


= 0 . 


(19) 


Due to the fact that is complex, it is split into a real- and imaginary part: 


ut := Re'I'(a;,t -|- At), vt := Im4'(a:,t -I- At), u := Re4'(a:,t), and v := Im\I'(a:,t). (20) 

The space and time discretised version of the Gross-Pitaevskii equation then becomes a system of two coupled 
non linear differential equations for the unknown functions Ut and vt- 


(ut-u - ^ (-A -F V) {vt +v)- ^ (ivt + u)^ -f {ut + u)^J (vt + v)\ 
yut - u -I- (-A -I- V) (ut -I- u) -I- ^ ^(ut -I- u)^ -I- {ut + u)^'j (ut + u)J 


In order to find the new wave function at time step t + At we use the standard Newton method (not to be 
confused with our constrained Newton method). Multiplying both components of (|?T|) with two different test 
functions and integrating over the volume leads to the weak formulation. Carrying out the the FEM discretisation 
we find for the right hand side 


where 


/ a-g, - f V6+Vg, -b fVb+g, - ^ ((a+)" + (6+)') b+g, \ 

+ ^Va+V/i, + ^Va+h, + ^ ((a+)" + (6+)") a+h,j ’ 


( 22 ) 


■=Utix^) + u^{x^), a :=u^{x^)-u^{x^),b+■=v^{x^) + v^{x^),b := v^{x^) - v{x^). (23) 


Here, again, we sum over all cells K and over all quadrature points x^ of each cell, is the weight of the 
quadrature formula evaluated at x^ times the Jacobian of the unit cell, gi are the test functions for the real part and 
hi for the imaginary part evaluated at position x^. The reader should be aware that the practical implementation 
is slightly different than (1^^ suggests. As a matter of fact, the form here has been chosen for better understanding. 
In practice, both discretised functions representing the imaginary and real part of the wave function are represented 
by one vector containing all degrees of freedom. More information about vector valued problems in the framework 
of deal.H can be found in the deal.H documentation. The sketch of the whole algorithm can be found in figure [3l 


3 Running the program 

3.1 Available programs 

The whole package consists of eight programs. Six are used for the simulations, three of which are responsible for 
computations of stationary states, and the remaining three are used for the real-time propagation. The last two 
programs serve for generating parameter files. 


7 








Figure 3: Flow chart of the real time propagation. Here we defined a+ := u^{x^) + u^{x^), a := u^{x^) — u^{x^), 
h+ := vKx^) + v\x^), b- := ) - v{x^). 


breed_1 

One-dimensional solver for the stationary Gross- 
Pitaevskii equation. 

breed 

Solver for the stationary Gross-Pitaevskii equation in 
Gartesian coordinates. 

breed_cs 

Same as above but for cylinder symmetry. 

rtprop-l 

One-dimensional solver for the real-time propagation. 

rtprop 

Solver for the real-time propagation of the wave function 
in Cartesian coordinates. 

rtprop_cs 

Same as above but for cylinder symmetry. 

gen.params 

Generates a directory with sub-directories containing pa¬ 
rameter files 

for different quantum numbers. 

gen_params_cs 

The same like above, just for cylindrical symmetry. 


For taking advantage of parallelisation, the 2D- and 3D simulations (breed, breed_cs, rtprop and rtprop.cs) should 
be called via mpirun -np [no of CPU cores] [program name]. Of course, single core execution is also possible. 
Note that in case of the ID simulations only single core operation is supported. When the atus-pro package is com¬ 
piled for the first time without having modified cmake options beforehand, the programs breed, breed_cs, rtprop 
and rtprop_cs operate in 2D-mode. Full 3D mode is supported by breed and rtprop. 









3.2 Setup of trapping potential and parameters 

The potential V{x) is specified inside the source code where pre-defined traps are: (i.) gravito-optical surface trap 
and (ii.) harmonic trap. The desired type of trap and the dimensionality of the problem has to be specified prior 
to compiling the package via CMake options (see subsection ibH) . This can be changed at any time and requires 
recompiling the code after having made the changes via the according CMake options. In practice, recompilation 
just takes a small amount of time. In order to use arbitrary traps, you can define them in the source code of the 
respective solver. For more details, it is suggested to take a closer look at the doxygen HTML documentation. Note 
that adjustments of all parameters related to physical or numerical quantities, e.g. the trap frequencies, do not 
require recompiling the code since this is done through the parameter file. 

3.3 Parameter file 

All programs expect a parameter file named "params. prm" where the values of the necessary parameters are set. 
It must be located in the folder from which the simulation runs are invoked and can be created via the programs 
themselves upon first invocation from within any directory. The parameter file is then automatically generated with 
default values and can be edited afterwards. Alternatively, you can use the programs gen_params or gen_params_cs 
which will generate a folder with sub folders labelled with respect to the quantum numbers of the initial guess 
functions. Each of those sub folders contains its own parameter file. This set of cases is meant to be used 
with breed and breed_cs. In case that the predefined set of values of the parameters should be changed, editing 
gen_params. cpp and gen_params_cs. cpp is required, followed by a recompilation of the atus-pro package. 


3.4 On the behaviour of the stationary states solvers: single vs. multiple solutions 

As already mentioned, the default strategy of the solvers breed_1, breed and breed_cs is to generate consecutive 
solutions w.r.t. increasing chemical potential fi = fig + This is the standard setting and cannot be changed via 
the parameter file. If one is interested in calculating only one stationary solution, then the source files breed_1 . cpp, 
breed.cpp or breed_cs.cpp need to be edited. Note that only the source file of the corresponding program, which 
you like to use for the calculation, has to be modified. In each of the mentioned files, in the main section of the code, 
the default statement solver. run2b() must be turned into a comment and solver. run() has to be uncommented. 
This is then followed by a recompilation of the atus-pro package. 


3.5 Output 


Stationary states When the programs breedJ, breed and breed_cs are invoked, the output on the terminal 
gives you various information about the status of the simulation. In the case of the test-run, which will be discussed 
in section IH the first step of the calculation generates an output like this : 


- 0 - 0000/ # - 0 - : number of iteration; - 0000/ : subdirectory where the calculation takes place 

Solving... 

Counter == 0 # number of iteration 

res == 4.69515e-07 # L2 norm of L2 gradient (ideally it would be zero), 

resp == 0.000293776 # difference of res: res(iteration before)-res(current) 

Ires / respI == 0.00159821 # can be used to estimate convergence rate 

mu == 5.7 # Chemical potential 

gs == 1 # Self interaction strength 

tl == 1.34114 # 1st factor for constrained Newton algorithm 

t2 == 0.997539 # 2nd factor for contrained Newton algorithm 

12 norm t == 1.67145 # L2-norm of vector t=(tl,t2) 

total no of cells == 154624 # total number of cells of all nodes 

total no of active cells == 115969 # total number of cells being used during integration 


Note that we added short explanations indicated by the comment symbol “#”. Concerning the output variable res, 
the L2-gradient is the left hand side of the stationary GPE in equation ([3]) , i.e. a good solution is characterised by 
a value close to zero. Here, the default value of the L2-gradient is set to 10“®. It can be adjusted in the params. prm 
file by means of the variable set epsilon. Upon reaching this value, the calculation for the actual mu is finished and 
the norm N of the wave function is printed to the terminal. Note that the self interaction strength gs corresponds 
to the value of 7 in eq. ©. 
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Dynamics Upon starting rtpropJ, rtprop and rtprop_cs the terminal outputs typically show the following 
information: 

min_cell_diameter = 0.0828641 
max_cell_dianneter = 10.6066 
dt/dx“2 == 1.45636 # Stability criterion 
t == 0 # Time-stamp of initial state 
N == 1.8432 # Number of particles 

p == 0, 0, 0 # Expectation value of momentum (p_x,p_y,p_z) 

pos == 2.74494, -3.61401e-16, 0 # Expectation value of position (x,y,z) 

t == 0 # Initial time of simulation 

m_res = 0.000265367 # L2-norm of L2-gradient of time-dependent GPE 
m_res = 4.78803e-09 


This case refers to the initial state in the directory 0000/ for the test-run in Cartesian coordinates. Note that the 
user can get a rough idea of the grid in use by means of the values min_cell_diameter and niax_cell_diameter. 
Concerning the stability criterion, which is based on the Courant condition, we would like to emphasise that it is 
generally recommended to have dt/dx^ < 1/2. This value can be influenced by adjusting the time step dt in the 
parameter file params. prm. The default value is set dt = le-2 and setting it to le-3 will lead in the presented 
case to dt/dx"^ = 0.145636. 

Note that the calculations for each time step are carried out until the residual m_res reaches 10“^®. This setting 
can only be adjusted in the program code. Using a smaller value for the residual than the predetermined one is 
not suggested since it would be below machine epsilon for double precision floating point numbers and leads the 
programs to perform calculations ad infinitum. 

Files All programs store the data in the deal.II native format and in the vtu format. The first format is used for 
storing the triangulation data and the wave function. For example, you can use a stationary solution from breed 
as input to the real- time propagation program rtprop. The second file format is meant to be used for graphical 
post-processing of the simulation data. We tested this extensively under ParaView (http://www.paraview.org/) 
and Visit (https://wci.llnl.gov/simulation/coniputer-codes/visit/). Moreover, both software packages are 
able to perform various data post processing tasks. 

Stationary states: In the case of stationary states the output is written into the files: (a) final.bin, 
final. bin. info and, e.g. final-00001 .vtu (two- and three-dimensions), whereas (b) in the one-dimensional case 
these are final.bin and, e.g. final-00001 .gnuplot. Note that in the latter case the graphical output format is 
best suited for Gnuplot (http: //www.gnuplot. info/). These file names are derived from the "filename" parameter 
entry in "params. prm". 

Moreover, a result file results.csv is written to the folder from which the simulation has been initiated. It 
contains the following data (only one entry is shown): 

Listing 1 : results.csv 

# tnu;gs;N;Counter;status 
5.7;1;1.8432;!;0 

Here, mu is the dimensionless chemical potential, gs is the nonlinearity 7, N is the norm of the resulting wave 
function, Counter refers to the number (Counter-f 1 ) of iterations which was needed for reaching convergence and 
finally status indicates whether the calculations have been successful (status = 0) or not (status = 1). 

Moreover, the status file log.csv will be created. In the case of single solution computations (see subsection 
lT4l) it is located in the folder from where the simulation has been started, whereas for multiple solutions, this file 
is found in each respective sub-directory. Typically, it has the following entries: 

Listing 2: log.csv 

# Counter;res;resp;I res / resp|;mu;gs;t1;t2;12 norm t;total no of cells;total no of active cells 
0;4.69514e-07;0.000293776;0.00159821;5.7;1;1.34114;0.997539;1.67145;154624;115969 

These are the same values like those discussed in subsection 13.51 

Real time propagation: For the dynamical case the following output files are generated (example: time step 
t = 0.1): (a) solution-0.100000. vtu, while (b) solution-0.100000.gnuplot, and so on. 
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Note that to initiate a run of the programs rtprop and rtprop.cs for two- and three-dimensional problems 
the files final.bin and final.bin.info are required. In contrast, for the time propagation in one dimension 
rtpropJ, no file final. bin. info is needed. In fact, it is not generated by breed_1. Note that for all real-time 
simulation programs the parameter file params. prm is needed. - In addition, rtprop and rtprop_cs create a file 
named solution.pvd which contains a temporal sequence of the results for intermediate time steps. This applies 
in case you like to produce a film of the real-time propagation with ParaView. 





Figure 4: Rescaled nonlinearity 7 = N ■ gs as function of the chemical potential fi. The solutions have been obtained in the 
testrun (section |4]) by the program breed for multiple solutions of the first excited state of a BEC in a GOST. Each fj,i is 
eigenvalue to the corresponding solution of the stationary Gross-Pitaevskii equation. 


3.6 Normalisation of the wave function and fixing the particle number 

As mentioned earlier in l2.21 the resulting wave functions from the stationary solvers are not normalised. In order to 
normalise them to one, the resulting norm of the wave function N can be used from the file results.csv together 
with the value for gs. That is, one selects for a particular chemical potential fXi the corresponding norm Ni. Here, 
the index i runs from 1 to the number of the last row in results.csv. Then the resulting nonlinearity can be 
set to — >■ 7 i = Ni • gs and the computed wave function for the pair (/ii,Ni) must be rescaled according to 

^ 4>i = In terms of the stationary GPE ([3]) the new quantities now fulfill the eigenvalue equation 

+ (25) 

with normalisation 

[ {cj)i{x))'^(f‘x = 1 . (26) 

Jv 

Now suppose that the initial (dimensionless) GPE ([3]) is the result of rescaling the (dimensionful) equation with 
the help of a typical length-scale L of the system, 

X ^ L ■ X, (j){x) ^ (j){x) ■ , (27) 


and the rescaled wave function is supposed to be normalised to one. In this case the dimensionless non-linearity 
term reads 


2mgs 


A-nh^as 

where Qs = - 

m 


(28) 


thus 


STra., 


7p = 


(29) 


Here as is the s-wave scattering length and m is the mass of the atom species. Because of the fact that the 
computations are made with non-normalised wave functions, 7 p cannot be identified with the initial nonlinearity 7 
of equation ©. However, by means of eq. ([25|l one can identify the following products: 


7 = 7 • N = 7p • Np, 


(30) 
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since the wave function (j){x) is normalised to one. Note that Np is the real, physical value of the particle number. 
By means of this one gets the particle number 


Np = —— 
^ SttOc 


L. 


(31) 


Alternatively, there exists the possibility to read the particle number Np directly from the table in results.csv. 
One demands 

7 = 7p, (32) 

which leads to N = Np and 7 = Sirag/L. Upon fixing the length-scale L, the value of 7 , that is gs, is determined 
and can be adjusted in the parameter file params. prm. From now on, the value N in results. csv has the meaning 
of a real particle number. 

If one wishes to perform simulations for a determined number of particles A®, it is suggested to compute first 
multiple solutions with different (see the testrun in section H. The resulting ( 7 i,/ii) pairs can be visualised by 
plotting the first column versus the product of the second and third column of the file results.csv. In this way 
one can get an idea about the functional relation between 7 ^ and Hi (see figure 2]). By means of equation (1311) 7 ° 
is determined by the desired particle number N^. Finally, the pair (/r'', 7 °) can be read (or interpolated) from the 
plot. 


3.7 Practical guidelines and hints 

Concerning the stationary solvers, if the external potential V{x) is too shallow, it could happen that the solver 
“jumps” between different solution branches. That means, the outcome of the solution becomes partly unpredictable, 
i.e. not controllable. Figuratively speaking, the wave function leaves the potential partially and gets delocalised, 
covering the whole simulation box. Thus the corresponding eigenfunctions are also possible solutions. In order to 
avoid this, one should scale the trap accordingly in order to get a deeper potential. 

Should the convergence behaviour be problematic, in many cases, it can be handled by modifying two parameters 
in the file params. prm. We suggest to decrease the increment A/r via the value of dmu or/and adjust the damping 
factor df of the contrained Newton algorithm, which is set per default to 1 and should be in the range df G [0,1]. 

As a general statement, try to avoid choosing an unsuitable spatial scale w.r.t. the typical grid scale. This could 
happen when the GPE is rescaled with a length scale much smaller than the average cell length. It is then highly 
probable that the solvers won‘t converge. In our simulations we usually take 10“® m as scaling length. Of course, 
this depends on the external potential in use. 



Figure 5: Three different solution of the stationary GPE with V{x,y) = 2;/2 -I- 2/^/4, 7 = 1 . (a) This is the ground state 
depicted together with the defanlt grid, (b) Results found by breed should look like this, if nothing is changed in the 
params. prm file, (c) An exemplary solution with a complicated structure. 
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(a) fi = 6.1, N « 225.25, 
QNl.x=0, QNl.y=0 
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(b) ti = 8.1, N « 265.76, 
QNl_x=1, QNl_y=1 



-1^ 


.i'dci 1 

(c) n= 13.1, iV« 378.14, 
QNl.x=4, QNl.y=3 


Figure 6 : Three different solution of the stationary GPE with V{x,y) = x^/i + i/^/4, 7 = 1 . a) This is the ground state 
depicted together with the default grid, (b) In case of BUILD_HTRAP=OFF a solution of a default run breed should look like 
this, (c) An exemplary solution with a complicated structure. 


4 Example of a test-run 

4.1 Stationary states 

A test-run may be performed within any directory by simply executing mpirun -np [no of CPU cores] breed. 
In case that the parameter file params.prm is missing, this will automatically generate one for a default scenario 
which is a gravito-optical surface trap in two dimensions (see equation ([1])). Then, before starting the calculations, 
breed creates ten subfolders which are sorted in ascending order with respect to the chemical potential /i. The 
value of each /r can be looked up in the file log.csv in each of the subfolders. Note that the simulation then runs 
immediately and generates stationary states in Cartesian coordinates for each of the ten predetermined values of /i. 
Here, the initial guess functions are chosen with quantum numbers {rix, Uy) = (1,1). Finally, the results are written 
into each of the mentioned subfolders into the files final.bin, final.bin. info and final-00001 .vtu (these are 
the default names). The latter serves for visualisation purposes and can be read, for example, by ParaView. As a 
visual example for this test-run, please see figures [5] and [6l 


4.2 Dynamics 

The previously generated stationary states can serve as initial states for the simulation of the time-evolution. For 
this, mpirun -np [no of CPU cores] rtprop can be initiated in any of the subfolders, where the files final.bin, 
final. bin. info and params. prm are present. In the default setup, the simulation runs until ttotai = 1 (dimensionless 
units) in steps of dt = le — 2 has been reached, where after each time step At = 0.1 an output file is written. This 
can be adjusted in the parameter file params. prm by means of the parameter NA accounting for the frequency of data 
output. In the case of this test-run ttotai — 1 and NA = 10 which translates to ttotai • NA = 10 output files in total. 
Another parameter which is adjusted in params.prm is NK specifies the number of intermediate steps between two 
file outputs. The results of the intermediate steps are printed in the terminal from which the simulation has been 
started. In the test-run the default value is NK = 10, therefore a total number of NA • NK = 100 steps is performed in 
the simulation, each with time step dt = le — 2. Thus, the total running time ttotai is determined by the product 
ttotai = NA • NK • dt. 


5 Conclusions and Outlook 

We have presented the C-F-F package atus-pro for computation of solutions of the stationary as well as the time- 
dependent Gross-Pitaevskii equation. The code incorporates finite element methods in order to speed up calculations 
for BECs with complex spatial structures. These could be invoked by excited states solutions, e.g. high oscillation 
quantum numbers of the guess function, as well as by elaborated three-dimensional external potential setups. 

Here we would like to point out briefly the main differences between our code and those described in m 
and [12]. For the solution of the non-stationary case we use the fully implicit Crank-Nicolson method, leading 
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to a system of non-linear equations. These are solved with a standard Newton method. In comparison, in the 
abovementioned works the semi-implicit Crank-Nicolson scheme is used and an operator splitting is carried out. 
This is the Alternating-Direction Implicit (ADI) method where the Laplace operator is split. Furthermore, this 
has been implemented for finite differences, whereas we use finite element methods. Finally, and this is the main 
difference, our stationary solver allows for calculation of excited states solutions. 

Excited topological modes have been also considered in [3S], where a damping gradient method is used. We like 
to point out in the following, how that method differs from our constrained Newton algorithm. Finding the ground 
state in the non linear case with gradient based methods (Quasi-Newton) should work always, but finding saddle 
point solutions (topological modes, non-ground state solutions) becomes difficult for high non linearities and high 
excitations. Finding the first or second excited state is possible with a good Quasi-Newton method with a good 
initial guess. The mathematical reason is that the Morse index increases for complicated solutions. The Morse 
index is the number of negative eigenvalues of the second variational derivative of equation (ED- This basically 
provides the information how many linearly independent descent directions are available. In the numerical world 
one can only deal with approximated solutions, which lie near critical points. Therefore, gradient based methods 
tend to miss critical points if the Morse index is too high. Sometimes they even simply miss it. The idea of the 
contrained Newton method goes back to Mountain pass or more general, to the Linking theorem (see the references 
in m)- The search is restricted to a special manifold and the saddle points lie on it. We cannot give a guarantee 
that all saddle points can be found, however we think that more solutions can be found than using other methods. 
This does not take into account numerical continuation methods. 

Numerical continuation is able to find more solutions, but there it is necessary to detect bifurcation points by 
evaluating, roughly speaking, the change in the Morse index. This is very time consuming for large problems in more 
then one spatial dimension. So at the end a even more powerful solver would incorporate a combination of both 
methods, that is including a constrained Newton method. A possible future extension is to implement numerical 
continuation. The aim of numerical continuation is to follow solution branches and detect possible branch points. 
Generally, a solution needs to be known beforehand. A Matlab tool box for numerical continuation which is based 
on this solution strategy is pde 2 path [26]. However, for our breed solvers it is not required to start with a known 
solution. In fact our code uses eigenfunctions of the corresponding linear problem just for convenience reasons. One 
could think of using the breed solvers to find more complicated solutions without knowing the branch. That is, by 
incorporating initial guess functions having a complex spatial structure. As suggested in m this can be utilised 
as a starting point for numerical continuation. In contrast to other procedures, it is then not necessary to detect 
branching points in order to find new branches. 

To enhance the range of physical applications of our code, an extension to dual-species BEG setups is intended. 
This is especially important in the frame of modelling experiments which test the Universality of Eree Fall (see, e.g. 
0 ). For this, an already existing two-species code will be adapted to the deal, ii library. 

Concerning technical extensions to the atus-pro code, it should be straightforward to extend for periodic 
boundary conditions. Furthermore, in order to be able to take into account a greater variety of domain geometries, 
our codes could be enhanced to incorporate the loading of grids, which can be generated by e.g. Gmsh or Open 
CASCADE. 

Finally, the implementation of adaptive grid refinement is envisaged. 
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A Installation of required libraries 

Since the atus-pro package depends on many extra libraries which have to be compiled, too, and this process 
sometimes can get cumbersome, we will give detailed compiling and installation instructions. 

The packages needed for a successful installation of our program code can be divided into two groups: the first 
contains packages which are included in the most common Linux distributions, namely 

1st Group 

• CMake (min 2.8) 

• Doxygen (optional) 

• Python (2.7), Bison, Flex 

• C, C++ 11 compliant compiler, Fortran 

• zlib (including also the development package) 

whereas in the second group software is listed, which - usually - has to be retrieved from external sources: 

2nd Group 

• MPl Implementation 

• GSL - GNU Scientific Library (1.16) 

• LAPAGK (3.5.0) 

• P4EST (1.1) 

• PETSG (3.6.1) 

• deal.ll (8.3.0). 

A.l General remarks on installation of dependencies 

As already mentioned, the packages of the 1st Group should be contained and be easily accessible in all common 
distributions of Linux. For those, no compiling is necessary, so we will only give installing instructions. On Ubuntu 
you can use the following command 

sudo apt-get install g++ gfortran cmake cmake-curses-gui bison flex doxygen. 

On Redhat (OpenSUSE) like Linux distributions you may use 

sudo yum (zypper) install gcc-c++ gcc-fortran doxygen cmake zlib-devel flex bison. 

In the following, we focus on the installing instructions for the packages listed in the 2nd Group. We assume 
that you possess root privileges in order to be able to install all dependencies in /opt. But of course you are free 
to install everything into your local home folder by adding ~ in front of /opt. 

For convenience, you should be able to copy and paste the commands from the boxes. In addition, the file 
INSTALL_OPTIONS.txt is included in the archive atus-pro_v 1 .O.tar.gz, where the install and compilation options 
for all the packages are also listed and can be copied to a terminal from which installation takes place. 

We suggest to extract all the downloaded packages listed in IA.2[ IA.3[ IA.4[ IA.5[ IA.6I and IA.7I to a temporary 
folder. In the following we assume that the folder name is temp and located in the top level of your home directory, 
i.e. the complete path would be then $H 0 ME/temp. 
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A.2 Installation of MPI 


You can skip this step if you have an MPI implementation already installed, which has also Fortran enabled. Here 
we present the installation of Open MPI, but other implementations of the MPI standard should work as well. 

It is important to check the lib folder naming of each library after each installation. Depending on your Linux 
distribution this may differ. It can be labelled as lib or lib 64 on 64-bit operating systems, whereas on 32-bit 
operating systems it should be always named as lib. You may check this via the command uname -a beforehand. 
Note that a wrong LD_LIBRARY_PATH will cause a failure of the compilation process. 

(0) cd $H0ME # start from top level of your home folder 

(1) mkdir temp # create temporary folder 

(2) cd temp 

(3) wget http://www.open-mpi.org/software/ompi/vl.8/downloads/openmpi-1 .8.5.tar.gz 

(4) tar xfvz openmpi-1.8.5.tar.gz 

(5) cd openmpi-1.8.5 

(6) ./configure —prefix=/opt/openmpi-1.8.5 —enable-mpi-fortran 

(7) make -j 

(8) sudo make install 

(9) export PATH=$PATH:/opt/openmpi-1.8.5/bin 

(10) export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/openmpi-l.8.5/lib64 

# alternativly: export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/openmpi-1.8.5/lib 

(11) check the path with: which mpicc 

(12) check the lib folder naming 


A.3 Installation of GSL 

The GNU Scientific Library [28], or short GSL, is required for the ansatz functions and for the search of the reference 
point in breed and breed_cs. 


(0) cd $H0ME/temp 

(1) wget ftp://ftp.gnu.org/gnu/gsl/gsl-l.16.tar.gz 

(2) tar xfvz gsl-1.16.tar.gz 

(3) cd gsl-1.16 

(4) ./configure —prefix=/opt/gsl-l.16 

(5) make CFLAGS="-march=native -03" # here we overload the default CFLAGS options 

(6) sudo make install 


A.4 Installation of LAPACK 

Lapack is needed by deal.ii for the single core direct solver called by the programs for the one-dimensional 
simulation breed_1 and rtprop_1 . Other available LAPACK versions may also work, but they have not been tested 
so far. 

(0) cd $H0ME/temp 

(1) wget http://www.netlib.Org/lapack/lapack-3.5.0.tgz 

(2) tar xfvz lapack-3.5.0.tgz 

(3) cd lapack-3.5.0 

(4) mkdir build 

(5) cd build 

(6) cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_Fortran_FLAGS_RELEASE="-march=native -fpic -03" \ 
-DCMAKE_INSTALL_PREFIX=/opt/lapack-3.5.0 .. 

(7) make -j 

(8) sudo make install 


A.5 Installation of p4est 

The p4est library requires MPI. The p4est library |2^ is responsible for the management of the parallel triangulation. 
It is important to mention that from now on the order of installation of the last three packages (including this one) 
is important and has to be performed in the order of the numbering of the subsections, i.e. install lA.Si IA.6I and 
finally IA.7I 
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(0) cd $H0ME/temp 

(1) wget http://p4est.github.io/release/p4est-1.1.tar.gz 

(2) tar xfvz p4est-1.1.tar.gz 

(3) cd p4est-1.1 

(4) ./configure —prefix=/opt/p4est-1.1 —enable-mpi —enable-shared —disable-vtk-binary —without-blas 

(5) make CFLAGS=”-03 -march=native'' 

(6) sudo make install 


A.6 Installation of PETSc 

The PETSc [30l [STJ [32] library requires MPI and a working internet connection. The PETSc library provides the 
deal.II library with routines which are required to solve systems of linear equations in parallel. We prefer the parallel 
direct solver MUMPS (http://mumps.enseeiht.fr). 

(0) cd $H0ME/temp 

(1) wget http://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.6.1 .tar.gz 

(2) tar xfvz petsc-3.6.1.tar.gz 

(3) cd petsc-3.6.1 

(4) export PETSC_ARCH=x86_64 # you can choose any arbitrary string you like 

(5) ./configure —prefix=/opt/petsc-3.6.1 —with-shared-libraries —with-x=0 —with-debugging=0 \ 

—with-mpi=1 —download-hypre=yes —download-fblaslapack=1 —download-scalapack —download-mumps \ 

—download-ptscotch 

(6) make PETSC_DIR=$H0ME/temp/petsc-3.6.1 PETSC_ARCH=x86_64 all 

(7) sudo make PETSC_DIR=$H0ME/temp/petsc-3.6.1 PETSC_ARCH=x86_64 install 


A.7 Installation of deal.II 

The deal.II library requires MPI and the preceding libraries. In order to speed up the compilation process you 
can optionally invoke ccmake . from within the folder $H0ME/temp/dealii-8.3.0/build (see the box below) and 
change the entry in CMAKE_BUILD_TYPE to Release. 

Moreover, if you want to benefit from an extra speed-up you can utilise more (or all) available CPU cores 
through executing make -j[number of CPU cores] in step (7). As a precaution, we do not recommend using just 
make -j (i.e. without explicitly determining the number of CPU cores) since the build process may slow down your 
computer extremely ! 

(0) cd $H0ME/temp 

(1) wget https://github.eom/dealii/dealii/releases/download/v8.3.0/dealii-8.3.0.tar.gz 

(2) tar xfvz dealii-8.3.0.tar.gz 

(3) cd dealii-8.3.0 

(4) mkdir build 

(5) cd build 

(6) cmake -DDEAL_II_WITH_UMFPACK=0N -DDEAL_II_WITH_LAPACK=0N -DLAPACK_DIR=/opt/lapack-3.5.0 \ 
-DPETSC_ARCH=x86_64 -DPETSC_DIR=/opt/petsc-3.6.1 -DP4EST_DIR=/opt/p4est-1.1 \ 

-DDEAL_II_WITH_THREADS=OFF -DDEAL_II_WITH_MPI=0N -DDEAL_II_WITH_HDF5=0FF \ 
-DCMAKE_INSTALL_PREFIX=/opt/deal.II-8.3.0 .. 

# after cmake has finished, it prints a summary of your configuration. The configuration should look 

# like the output in the box below. 

(7) make # alternatively try make -j[no of epu cores] 

(8) sudo make install 

After step (7), upon successful compilation, the terminal output should look like this for the user “johndoe”: 

### 

# 

# deal.II configuration: 

# CMAKE_BUILD_TYPE: 

# BUILD_SHARED_LIBS: 

# CMAKE_INSTALL_PREFIX: 

# CMAKE_SOURCE_DIR: 

# 

# CMAKE_BINARY_DIR: 


DebugRelease 

ON 

/opt/deal.II-8.3.0 

/home/johndoe/temp/dealii-8.3.0 

(version 8.3.0) 

/home/johndoe/temp/dealii-8.3.0/build 
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GNU 4.8.3 on platform Linux i686 
/usr/bin/c++ 


# CMAKE_CXX_COMPILER: 

# 

# 

# Configured Features (DEAL_II_ALLOW_BUNDLED = ON, DEAL_II_ALLOW_AUTODETECTION = ON); 

# ( DEAL_II_WITH_64BIT_INDICES = OFF ) 

# ( DEAL_II_WITH_ARPACK = OFF ) 

# DEAL_II_WITH_BOOST set up with external dependencies 

# ( DEAL_II_WITH_BZIP2 = OFF ) 

# DEAL_II_WITH_CXX11 = ON 

# ( DEAL_II_WITH_CXX14 = OFF ) 

# ( DEAL_II_WITH_HDF5 = OFF ) 

# DEAL_II_WITH_LAPACK set up with external dependencies 

# ( DEAL_II_WITH_METIS = OFF ) 

# DEAL_II_WITH_MPI set up with external dependencies 

# DEAL_II_WITH_MUPARSER set up with bundled packages 

# ( DEAL_II_WITH_NETCDF = OFF ) 

# ( DEAL_II_WITH_OPENCASCADE = OFF ) 

# DEAL_II_WITH_P4EST set up with external dependencies 

# DEAL_II_WITH_PETSC set up with external dependencies 

# ( DEAL_II_WITH_SLEPC = OFF ) 

# ( DEAL_II_WITH_THREADS = OFF ) 

# ( DEAL_II_WITH_TRILINOS = OFF ) 

# DEAL_II_WITH_UMFPACK set up with bundled packages 

# DEAL_II_WITH_ZLIB set up with external dependencies 

# 

# Component configuration: 

# ( DEAL_II_COMPONENT_DOCUMENTATION = OFF ) 

# DEAL_II_COMPONENT_EXAMPLES 

# DEAL_II_COMPONENT_MESH_CONVERTER 

# ( DEAL_II_COMPONENT_PACKAGE = OFF ) 

# ( DEAL_II_COMPONENT_PARAMETER_GUI = OFF ) 

# 

# Detailed information (compiler flags, feature configuration) can be found in detailed.log 

# 

# Run $ make info to print a help message with a list of top level targets 

# 

### 

— Configuring done 

— Generating done 

— Build files have been written to: /home/johndoe/temp/dealii-8.3.0/build 
johndoeOlinux:'/temp/dealii-8.3.0/build> 


B Building the atus-pro package 

Before we can build the atus-pro package from the source it is necessary to setup the environmental variables. The 
PATH variable should point to the folder where the MPI binaries are located. In order to be capable to run our 
programs, in addition, it must include $H 0 ME/bin. The LD_LIBRARY_PATH variable should contain all paths pointing 
to all library locations. Our CMake configuration uses this path to search for the library files. 

We suggest to create a file (e.g. .my_bashrc) located in the home folder $H 0 ME with the contents given in 
the box below. This file can then be used to set up the paths by invoking the source command, i.e. source 
$H 0 ME/. myjDashrc. In general, we recommend the use of environment modules (http : //modules. sourceforge . net/) 
which is common on distributed memory machines. 

Listing 3: .my_bashrc 

export PATH=$PATH:$HOME/bin:/opt/openmpi-1 .8.5/bin 

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/openmpi-1.8.5/lib64:/opt/gsl-1.16/lib64:/opt/p4est-l.1/lib64: \ 
/opt/petsc-3.6.1/lib:/opt/deal.11-8.3.0/lib 

# alternatively 

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/openmpi-1.8.5/lib:/opt/gsl-1.16/lib:/opt/p4est-1.1/lib: \ 
/opt/petsc-3.6.1/lib:/opt/deal.11-8.3.0/lib 
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At this point we remind you to check the lib folder naming above in LD_LIBRARY_PATH. On 64-bit operating system 
it may be lib or lib 64 , whereas on 32-bit operating systems it should be always named lib ! For this, it is 
recommended to check the real installation location of the libraries listed in LD_LIBRARY_PATH. Note that by using 
the command uname -a you can verify if you are using a 32- or 64-bit operating system. 

Before starting to compile atus-pro, we assume that you have already copied the archive atus-pro_v1 .0. tar. gz 
into your $H 0 ME/temp folder ! 

(0) source $H0ME/.my_bashrc # establish the paths of the previously compiled libraries 

(1) cd $H0ME/temp 

(2) tar xfvz atus-pro_vl.0.tar.gz 

(3) cd atus-pro_vl.0 

(4) mkdir build 

(5) cd build 

(6) cmake . . 

(7) make -j # after the build process you can find the binaries in $H0ME/bin 

(8) make doc # optionally: build the doxygen documentation 

Within the build folder you can view all CMake options and preferences through invoking the ccmake . com¬ 
mand. If CMake has doxygen detected, then you can build the HTML documentation via the command make doc. 
The documentation is then accessible through the index.html located in the HTML folder. This documentation 
contains a more detailed description of the algorithms. 

B.l Important CMake options 

The CMake options can be accessed by executing ccmake . in the folder $H 0 ME/temp/atus-pro_v 1 . 0 /build. Editing 
the following listed options is necessary in case you decide to change the default settings (i.e. switch from 2D to 
3D problems). In order for the changes to take effect, this must be done prior to compiling the atus-pro package. 

If you build the package for the first time without changing these options, then the initial setup of the external 
potential is a two-dimensional gravito-optical surface trap in accordance with equations o and ©■ One can switch 
to a purely harmonic trap by setting BUILD_HTRAP ON. The values of each trapping frequency (wa,, cuy, Wz) are then 
fixed in the parameter file params.prm. Of course, if three-dimensional simulations are to be made, BUILD_3D must 
be set to ON before compiling atus-pro. 


CMake option 

default 

value 


BUILD_ 3 D 

OFF 

If this option is set to ON, then our 
C-F-F templates are instantiated for 
three spatial dimensions. 

BUILDDOCUMENTATION 

AUTO 

If Latex and doxygen are available on 
your computer, then the build target 
for the documentation is automati¬ 
cally enabled . 

BUILD_HTRAP 

OFF 

If this option is set to ON, then a 
harmonic trap is switched on. 

BUILDNEHARI 

ON 

Influences the initial point in the 
function space in breed and breed_cs. 
For a detailed explanation look into 
the doxygen documentation. 

BUILD_VARIANT 2 

ON 

Influences the search of the reference 
point in the function space in breed 
and breed_cs. For a detailed expla¬ 
nation look into doxygen HTML doc¬ 
umentation. 
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B.2 Overview of parameters of the parameter file params. prm 


NA 

Frequency of data output 

NK 

Number of intermediate steps 

Ndmu 

Number of A/r steps 

dmu 

A/i 

df 

Damping factor for the Newton method 

dt 

Time step At 

epsilon 

Termination threshold for stationary solvers. 

filename 

File name of the wave function for the initial time step 

guess_fct 

Set (j)Q manually (see doxygen documentation) 

ti 

Initial value for point in function space. If BUILDJvlEHARI = OFF, ^ ti and ^ ti 

global_refinements 

Level of global mesh refinements for the inital mesh 

xMax 

Max. value of first coordinate of simulation box 

xMin 

Min. value of first coordinate of simulation box 

yMax 

Max. value of second coordinate of simulation box 

yMin 

Min. value of second coordinate of simulation box 

zMax 

Max. value of third coordinate of simulation box 

zMin 

Min. value of third coordinate of simulation box 

QN 1 _x 

Quantum number corresponding to the first coordinate 

QNI.y 

Quantum number corresponding to the second coordinate 

QN 1 _z 

Quantum number corresponding to the third coordinate 

gf 

Gravitational acceleration 

gs 

Self interaction parameter 7 

omega_x 

UJ^ 

omega.y 

COy or CJr 

omega_z 

iOz 

t 

Starting time for the real time propagation 
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